library(fixest) # to demean variables
library(multiwayvcov) # to cluster SEs
library(sandwich) # to get HC2 SEs
library(stargazer) # to generate tables
library(mvtnorm) # to simulate betas and vcov matrix
library(car)# to transform DV

rm(list = ls())

# load data
load(file = 'chData.RData')

# Put each dataset in a different object
ch_reform <- ch_list[[1]]
ch_refTerm <- ch_list[[2]]
ch_pre_post <- ch_list[[3]]
ch_refTermMChanges <- ch_list[[4]]
ch_refTermParty <- ch_list[[5]]
ch_refTermMChangesParty <- ch_list[[6]]

# Exclude "Gonzalo Fuenzalida" and 'Javier Macaya' because their district changed after the reform
ch_reform <- subset(ch_reform, !(bill_author %in% c("Gonzalo Fuenzalida", 'Javier Macaya')))
ch_refTerm <- subset(ch_refTerm, !(bill_author %in% c("Gonzalo Fuenzalida", 'Javier Macaya')))
ch_pre_post <- subset(ch_pre_post, !(bill_author %in% c("Gonzalo Fuenzalida", 'Javier Macaya')))
ch_refTermMChanges <- subset(ch_refTermMChanges, !(bill_author %in% c("Gonzalo Fuenzalida", 'Javier Macaya')))
ch_refTermParty <- subset(ch_refTermParty, !(bill_author %in% c("Gonzalo Fuenzalida", 'Javier Macaya')))
ch_refTermMChangesParty <- subset(ch_refTermMChangesParty, !(bill_author %in% c("Gonzalo Fuenzalida", 'Javier Macaya')))


# Generate DVs
ch_refTerm$women_pct_logit <- logit(ch_refTerm$women_pct/100)
ch_refTermMChanges$women_pct_logit <- logit(ch_refTermMChanges$women_pct/100)
ch_reform$women_pct_logit <- logit(ch_reform$women_pct/100)
ch_pre_post$women_pct_logit <- logit(ch_pre_post$women_pct/100)

####### Models using M ##########
ch_refTerm$magAnalysisWoman <- ch_refTerm$magAnalysis * ch_refTerm$female
ch_reform$magAnalysisWoman <- ch_reform$magAnalysis * ch_reform$female
ch_pre_post$magAnalysisWoman <- ch_pre_post$magAnalysis * ch_pre_post$female
ch_refTermMChanges$magAnalysisWoman <- ch_refTermMChanges$magAnalysis * ch_refTermMChanges$female

# M1: G*M C1,C2,C3 (all data)
data_m1 <- fixest::demean(women_pct_logit + magAnalysis + female + magAnalysisWoman
	~ session, data = ch_refTerm)
data_m1 <- cbind(data_m1, bill_author = ch_refTerm$bill_author)
m1 <- lm(women_pct_logit ~ magAnalysis + female + magAnalysisWoman - 1, data = data_m1)
vcov1 <- cluster.vcov(m1, cluster = data_m1$bill_author)
s1  <- sqrt(diag(vcov1))

# M2: G*M C1,C2,C3 (M changes at the middle of C2)
data_m2 <- fixest::demean(women_pct_logit + magAnalysis + female + magAnalysisWoman
	~ session + reform, data = ch_refTermMChanges)
data_m2 <- cbind(data_m2, bill_author = ch_refTermMChanges$bill_author)
m2 <- lm(women_pct_logit ~ magAnalysis + female + magAnalysisWoman - 1, data = data_m2)
vcov2 <- cluster.vcov(m2, cluster = data_m2$bill_author)
s2  <- sqrt(diag(vcov2))


# M3: G*M C1,C3
data_m3 <- fixest::demean(women_pct_logit + magAnalysis + female + magAnalysisWoman
	~ session, data = subset(ch_refTerm, session %in% c('2010-2014', '2018-2022')))
data_m3 <- cbind(data_m3, bill_author = subset(ch_refTerm, session %in% c('2010-2014', '2018-2022'))$bill_author)
m3 <- lm(women_pct_logit ~ magAnalysis + female + magAnalysisWoman - 1, data = data_m3)
vcov3 <- cluster.vcov(m3, cluster = data_m3$bill_author)
s3  <- sqrt(diag(vcov3))

# M4: C2
inSample <- names(table(ch_pre_post$bill_author)[table(ch_pre_post$bill_author) == 2])
ch_pre_post <- subset(ch_pre_post, bill_author %in% inSample)
data_m4 <- fixest::demean(women_pct_logit + magAnalysis + female + magAnalysisWoman
	~ bill_author, data = ch_pre_post)
data_m4 <- cbind(data_m4, bill_author = ch_pre_post$bill_author)
m4 <- lm(women_pct_logit ~ magAnalysis + magAnalysisWoman - 1, data = data_m4)
vcov4 <- cluster.vcov(m4, cluster = data_m4$bill_author)
s4  <- sqrt(diag(vcov4))

# M5: G*M C1,C2,C3 (only repeated legislators)
data_m5 <- fixest::demean(women_pct_logit + magAnalysis + magAnalysisWoman
	~ bill_author + session, data = subset(ch_refTerm, allSessions == 1))
data_m5 <- cbind(data_m5, bill_author = subset(ch_refTerm, allSessions == 1)$bill_author)
m5 <- lm(women_pct_logit ~ magAnalysis + magAnalysisWoman - 1, data = data_m5)
vcov5 <- cluster.vcov(m5, cluster = data_m5$bill_author)
s5  <- sqrt(diag(vcov5))

# M6: G*M C1, C2, C3 (M changes at the middle of C2) (only repeated legislators)
data_m6 <- fixest::demean(women_pct_logit + magAnalysis + magAnalysisWoman
	~ session + reform + bill_author, data = subset(ch_refTermMChanges, allSessions == 1))
data_m6 <- cbind(data_m6, bill_author = subset(ch_refTermMChanges, allSessions == 1)$bill_author)
m6 <- lm(women_pct_logit ~ magAnalysis + magAnalysisWoman - 1, data = data_m6)
vcov6 <- cluster.vcov(m6, cluster = data_m6$bill_author)
s6  <- sqrt(diag(vcov6))

# M7: G*M C1,C3 (only repeated legislators)
data_m7 <- fixest::demean(women_pct_logit + magAnalysis + magAnalysisWoman
	~ bill_author + session, data = subset(ch_refTerm, session %in% c('2010-2014', '2018-2022') & l10_l18 == 1))
data_m7 <- cbind(data_m7, bill_author = subset(ch_refTerm, session %in% c('2010-2014', '2018-2022') & l10_l18 == 1)$bill_author)
m7 <- lm(women_pct_logit ~ magAnalysis + magAnalysisWoman - 1, data = data_m7)
vcov7 <- cluster.vcov(m7, cluster = data_m7$bill_author)
s7  <- sqrt(diag(vcov7))

# M8: C3
data_m8 <- subset(ch_refTerm, session == '2018-2022')
m8 <- lm(women_pct_logit ~ magAnalysis + female + magAnalysisWoman, data = data_m8)
vcov8 <- vcovHC(m8, 'HC2')
s8  <- sqrt(diag(vcov8))



# Table F.25
stargazer(m1, m2, m3, m4, m5, m6, m7, m8, 
	se = list(s1, s2, s3, s4, s5, s6, s7, s8), 
	covariate.labels = c('M', 'Woman', 'M x Woman'), 
	add.lines = list(c('FE by Legislative Term', 'Yes', 'Yes', 'Yes', 'No', 'Yes', 'Yes', 'Yes', 'No'),
		c('FE by Reform', 'No', 'Yes', 'No', 'No', 'No', 'Yes', 'No', 'No'),
		c('FE by Legislator', 'No', 'No', 'No', 'Yes', 'Yes', 'Yes', 'Yes', 'No')),
	dep.var.labels = "Bills' Portfolio on Women's Issues (%)",
	omit.stat = c("ser",'f'),
	star.cutoffs = c(0.10, 0.05, 0.01),
	no.space = TRUE,
	single.row = FALSE,	
	notes.append = FALSE,
	title = "Association Between Legislative Portfolio, Gender, and District Magnitude--2014-2021  Cámara de Diputadas y Diputados de Chile",
	label = 't:mag'
)

#################### Substantive Effect ########################
## Scenarios: Women, M = 2
W2 <- data.frame(mag = 2, female = 1, magAnalysisWoman = 2 * 1)
## Scenarios: Women, M = 3 (post-reform)
W3 <- data.frame(mag = 3, female = 1, magAnalysisWoman = 3 * 1)
## Scenarios: Men, M = 2
M2 <- data.frame(mag = 2, female = 0, magAnalysisWoman = 2 * 0)
## Scenarios: Men, M = 3 (post-reform)
M3 <- data.frame(mag = 3, female = 0, magAnalysisWoman = 3 * 0)
# Women and Men Data
Wdata <- rbind(W2, W3)
Mdata <- rbind(M2, M3)

## Simulate betas and vcov
set.seed(123)
sims1 <- rmvnorm(n = 1000, mean = coef(m1), sigma = vcov1)
sims2 <- rmvnorm(n = 1000, mean = coef(m2), sigma = vcov2)
sims3 <- rmvnorm(n = 1000, mean = coef(m3), sigma = vcov3)
sims4 <- rmvnorm(n = 1000, mean = coef(m4), sigma = vcov4)
sims5 <- rmvnorm(n = 1000, mean = coef(m5), sigma = vcov5)
sims6 <- rmvnorm(n = 1000, mean = coef(m6), sigma = vcov6)
sims7 <- rmvnorm(n = 1000, mean = coef(m7), sigma = vcov7)
sims8 <- rmvnorm(n = 1000, mean = coef(m8), sigma = vcov8)

## Calculate Predicted Values
# Women
pred1_W <- apply(plogis(sims1 %*% t(Wdata[2,])) - plogis(sims1 %*% t(Wdata[1,])), 2, quantile, c(0.975, 0.5, 0.025)) *100
pred2_W <- apply(plogis(sims2 %*% t(Wdata[2,])) - plogis(sims2 %*% t(Wdata[1,])), 2, quantile, c(0.975, 0.5, 0.025)) *100
pred3_W <- apply(plogis(sims3 %*% t(Wdata[2,])) - plogis(sims3 %*% t(Wdata[1,])), 2, quantile, c(0.975, 0.5, 0.025)) *100
pred4_W <- apply(plogis(sims4 %*% t(Wdata[2, -2])) - plogis(sims4 %*% t(Wdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025)) *100
pred5_W <- apply(plogis(sims5 %*% t(Wdata[2, -2])) - plogis(sims5 %*% t(Wdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025)) *100
pred6_W <- apply(plogis(sims6 %*% t(Wdata[2, -2])) - plogis(sims6 %*% t(Wdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025)) *100
pred7_W <- apply(plogis(sims7 %*% t(Wdata[2, -2])) - plogis(sims7 %*% t(Wdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025)) *100
pred8_W <- apply(plogis(sims8 %*% t(cbind(1, Wdata)[2,])) - plogis(sims8 %*% t(cbind(1, Wdata)[1,])), 2, quantile, c(0.95, 0.5, 0.025)) *100

# Men
pred1_M <- apply(plogis(sims1 %*% t(Mdata[2,])) - plogis(sims1 %*% t(Mdata[1,])), 2, quantile, c(0.975, 0.5, 0.025)) *100
pred2_M <- apply(plogis(sims2 %*% t(Mdata[2,])) -plogis(sims2 %*% t(Mdata[1,])), 2, quantile, c(0.975, 0.5, 0.025)) *100
pred3_M <- apply(plogis(sims3 %*% t(Mdata[2,])) - plogis(sims3 %*% t(Mdata[1,])), 2, quantile, c(0.975, 0.5, 0.025)) *100
pred4_M <- apply(plogis(sims4 %*% t(Mdata[2, -2])) - plogis(sims4 %*% t(Mdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025)) *100
pred5_M <- apply(plogis(sims5 %*% t(Mdata[2, -2])) - plogis(sims5 %*% t(Mdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025)) *100
pred6_M <- apply(plogis(sims6 %*% t(Mdata[2, -2])) - plogis(sims6 %*% t(Mdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025)) *100
pred7_M <- apply(plogis(sims7 %*% t(Mdata[2, -2])) - plogis(sims7 %*% t(Mdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025)) *100
pred8_M <- apply(plogis(sims8 %*% t(cbind(1, Mdata)[2,])) - plogis(sims8 %*% t(cbind(1, Mdata)[1,])), 2, quantile, c(0.975, 0.5, 0.025)) *100


# Building blocks of table F26
bill_hat_women <- c(round(pred1_W[2,], 2), round(pred2_W[2,], 2), round(pred3_W[2,], 2),
	round(pred4_W[2,], 2), round(pred5_W[2,], 2), round(pred6_W[2,], 2),
	round(pred7_W[2,], 2), round(pred8_W[2,], 2))
ci_women <- c(paste0('(', round(pred1_W[3,], 2), ', ', round(pred1_W[1,], 2), ')'), 
	paste0('(', round(pred2_W[3,], 2), ', ', round(pred2_W[1,], 2), ')'), 
	paste0('(', round(pred3_W[3,], 2), ', ', round(pred3_W[1,], 2), ')'), 
	paste0('(', round(pred4_W[3,], 2), ', ', round(pred4_W[1,], 2), ')'), 
	paste0('(', round(pred5_W[3,], 2), ', ', round(pred5_W[1,], 2), ')'), 
	paste0('(', round(pred6_W[3,], 2), ', ', round(pred6_W[1,], 2), ')'), 
	paste0('(', round(pred7_W[3,], 2), ', ', round(pred7_W[1,], 2), ')'), 
	paste0('(', round(pred8_W[3,], 2), ', ', round(pred8_W[1,], 2), ')'))


bill_hat_men <- c(round(pred1_M[2,], 2), round(pred2_M[2,], 2), round(pred3_M[2,], 2),
	round(pred4_M[2,], 2), round(pred5_M[2,], 2), round(pred6_M[2,], 2),
	round(pred7_M[2,], 2), round(pred8_M[2,], 2))
ci_men <- c(paste0('(', round(pred1_M[3,], 2), ', ', round(pred1_M[1,], 2), ')'), 
	paste0('(', round(pred2_M[3,], 2), ', ', round(pred2_M[1,], 2), ')'), 
	paste0('(', round(pred3_M[3,], 2), ', ', round(pred3_M[1,], 2), ')'), 
	paste0('(', round(pred4_M[3,], 2), ', ', round(pred4_M[1,], 2), ')'), 
	paste0('(', round(pred5_M[3,], 2), ', ', round(pred5_M[1,], 2), ')'), 
	paste0('(', round(pred6_M[3,], 2), ', ', round(pred6_M[1,], 2), ')'), 
	paste0('(', round(pred7_M[3,], 2), ', ', round(pred7_M[1,], 2), ')'), 
	paste0('(', round(pred8_M[3,], 2), ', ', round(pred8_M[1,], 2), ')'))

tableF26 <- rbind("Women "= bill_hat_women, " " = ci_women, 
	"Men "= bill_hat_men, " " = ci_men)
colnames(tableF26) <- paste0('Model ', 1:8)
stargazer(tableF26,
	title = "Predicted percentage point change in the cosponsorship portfolio on women’s issues when M increases by one")

